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ABSTRACT 

We analyze the orbits of stars and dark matter particles in the halo of a disk galaxy formed in a cosmological 
hydrodynamical simulation. The halo is oblate within the inner ~ 20 kpc and triaxial beyond this radius. About 
43% of orbits are short axis tubes - the rest belong to orbit families that characterize triaxial potentials (boxes, 
long-axis tubes and chaotic orbits), but their shapes are close to axisymmetric. We find no evidence that the self- 
consistent distribution function of the nearly oblate inner halo is comprised primarily of axisymmetric short- 
axis tube orbits. Orbits of all families, and both types of particles are highly eccentric with mean eccentricity 
> 0.6. We find that randomly selected samples of halo stars show no substructure in "integrals of motion" 
space. However individual accretion events can be clearly identified in plots of metallicity versus formation 
time. Dynamically young tidal debris is found primarily on a single type of orbit. However, stars associated 
with older satellites become chaotically mixed during the formation process (possibly due to scattering by 
the central bulge and disk, and baryonic processes), and appear on all four types of orbits. We find that 
the tidal debris in cosmological hydrodynamical simulations experiences significantly more chaotic evolution 
than in collisionless simulations, making it much harder to identify individual progenitors using phase space 
coordinates alone. However by combining information on stellar ages and chemical abundances with the orbital 
properties of halo stars in the underlying self-consistent potential, the identification of progenitors is likely to 
be possible. 

Subject headings: (Cosmology): dark matter, Galaxy: evolution, Galaxy: formation, Galaxy: halo, Galaxy: 
kinematics and dynamics, Methods: numerical 



1. INTRODUCTION 



In the picture proposed by |Searle & Zinnl ( |1978| ) the stel- 
lar halo of the Milky Way (hereafter MW) formed from the 
merger of protogalactic fragments consisting of dark matter 
(hereafter DM), gas, and the first generations of stars which 
formed in high density peaks in the early Universe. An ex- 
tended period of "late infall" and tidal shredding of satellite 
galaxies has continued to build up the halo even at the present 
time. In the past decade, resolved-star surveys have provided 
strong observational evidence in supp ort of this picture. Nu- 
merous coherent stellar tidal streams ([Newberg et al. 2002; 
|Yanny et aL||2003[ [Belokurov et al. 12006) and tens of ultra- 
faint dwarf spheroidal galaxies ( Willman et al. 2005 ; Zucker 
|et alJ2006{ |Belokurov et al.|2007| ) - possibly progenitors of 
the tidal streams - have been discovered in large surveys. 
Their discovery is broadly consistent with the predictions of 
galaxy formation simul ations in the cold dark matter (here- 



after ACDM) paradigm (Bu llock & Johnst on 2005 ; Bell et al. 
|2008) . 

In this paper we focus on the orbital properties of halo stars 
and dark matter particles in a disk galaxy produced in an 
A/-body+SPH simulation run in a fully cosmological context 
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( Stin son et al.||2010| part of the McMaster Unbiased Galaxy 
Simulations, hereafter MUGS). The MUGS sample consists 
of 16 galaxies (both ellipticals and spirals) with host halo 
masses of ~ 10 12 M selected in an unbiased way and resim- 
ulated with volume renormalization techniques. The MUGS 
simulation has higher mass resolution and smaller spatial 
gravitational softening than recent simulations which have ex- 
amined similar issues, but it still suffers from the over cool- 
ing problem that has plagued most cosmological simulations 
(however th i s situation has recently improved, e.g. [G o vernato 
eF^]j20TQl jGuedes et alT[|20TT| |Zemp et al.||2012[ |Stinson 
et al.|20l2p . The subject or our investigation is the most disk- 



dominated massive spiral galaxy in the MUGS sample, which 
happens to have a mass, and radial scale length similar to that 
of the Milky Way. 

We explore two questions in this paper (a) what is the shape 
of the dark matter halo of this disk galaxy as a function of ra- 
dius and what types of orbits self-consistently produce this 
shape? (b) is it possible to use phase space coordinates to 
associate halo stars with individual progenitor satellites in a 
fully cosmological context where baryonic physics may have 
altered the dynamics non-adiabatically? These are important 
questions in the context of modeling the Milky Way with stel- 
lar orbits and will become increasing important in the next 
few years with the advent of numerous resolved star surveys. 
One of the objectives of current and future surveys is to de- 
termine the global shape and density profile of the Milky Way 
halo either by modelin g the kinematics of local field halo stars 
( Loebman et al. 12012] ), and/ or combining this with models of 
individual tidal streams (e.g. Johnston et al.|1999"t|Law & Ma- ' 
jewski|2010| ). 



The motivation for our investigation of the first question is 
that although cosmological Af-body simulations produce dark 
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matter halos that are triaxial or prolate ( Jing & Suto 2002]), it 
has been known for over a decade that the shapes of DM halos 
are altered when baryons cool and condense at their centers. 
With baryons, halos become spherical or axisymmetric within 
the inner one third o f the virial radius, but remain triaxial or 
prolate at larger rad ii ( Dubinski & Carlberg 1991 ; Kazantzidis 
eTaL]|2004l IDebattista et al ||2008| ). If a gravitational poten- 



tial is exactly oblate and axisymmetric, all orbits in the sys 
tern are axisymmetric short axis tubes (BT08). The prevail 
ing view is that as a triaxial system becomes more oblate due 
to chaotic mixing or other processes, an increasing fraction 
of the box orbits and long axis tube orbits that constitute the 
"back bone" of the tria xial system are co nverted to axisym- 
metric short-axis tubes (Me rritt & Va lluri 1996 ). Although it 
is generally assumed that the self-consistent distribution func- 
tion of a nearly oblate system may be assumed to be com- 
prised exclusively of axisymmetric short-axis tubes, it has 
been shown that even small deviations from axisymmetry can 
result in a large fra ction of box orbits in the in ner regions of 
the potential ( [van d en Bosc h & de Ze euw 2009). Using con- 
trolled simulations in which a rigid particle disk grows adia- 



batically inside a triaxial halo , | Valluri et al. 1( 2010) (hereafter 
|VT0| ) and |Valluri eTaL] ( [20T2] ) (hereafter |V12j ) showed that al- 
though the halos become more oblate in the inner regions, a 
significant fraction of the orbit population retains its original 
orbital characteristics (e.g. box orbits remain box orbits and 
long axis tubes remain long-axis tubes). However, as the po- 
tential changes shape adiabatically, the halo orbits with small 
pericenter radii also become "rounder", but a large fraction 
do not change their orbital type because their orbital integrals 
of motion (or rather orbital actions) are adiabatic invariants. 
When a disk galaxy is formed in more realistic hierarchical 
galaxy formation simulations, the evolut ion is no longer adi- 
abatic. Although several recent studies ([Tissera et al.pOlO 



; al.|2010j|Zemp et al,|2012> find that the shapes 
r and stellar halos in cosmological simulations 
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of dark matter I 
are also nearly oblate, effects of baryonic processes on the or- 
bits of halo stars and dark matter particles in Milky Way sized 
disks in this co ntext is less well understood. A recent study 
( Bry an et al.|2012| ) suggests that orbital properties of stars and 
dark matter particles in 1O 13 M halos may be somewhat sen- 
sitive to the feedback prescriptions employed. Knowing how 
the shape of a cosmological halo varies with radius and how 
DM and stellar orbits respond can provide vital insights into 
the shape and formation history of our own halo. 

The motivation for investigating the second question is to 
study the impact of baryonic processes on "Galactic Archeol- 
ogy" with halo stars. Although halo stars probably constitute 
no more than 1% of the stellar mass of the Galaxy, current 
and future surveys are focused on uncovering their kinemati- 
cal, spatial and abundance distributions. This is because it has 
been argued that the time it takes for the integrals of motion of 
the orbits of halo stars to change may be longer than the age 
of the Galaxy, allowing correlations between the abundances 
and kinematics of halo stars to serve as a "fossil" record of the 
formation history of the Galaxy. It has been shown that de- 
constructing the accretion history of the halo (i.e. identifying 
individual accretion events) is possib le using orbital integrals 
of motion in controlled simulatio ns ( Helm i & de Zee uw 2000; 



McM illan & Binne y 2008; Gom ez et aLpO lO). However, the 
effects on the orbital integrals of motion of hierarchical evolu- 
tion accompanied by dissipation, star formation and feedback 
from baryons, have yet to be examined. 



In addition, there has been significant recent effort devoted 
to interpreting the kinematic, metallicity and orbital distribu- 
tions of halo stars to infer the formation history of the stellar 
halo. Theoretical studies have identified three possible ori- 
gins for stars that currently reside in the halo: (a) they formed 
in a satellite galaxy outside the virial radius and were sub- 



sequently accreted (the so called "accreted" halo, Bullock & 
Johnston 2005 ), (b) they formed inside a satellite after it is 
accrete d by the main galax y {"in situ halo stars" or "endo- 
debris" Tissera et al. 2013 ), (c) they formed in the disk and 
were then kicked in to the halo by infalling satellites ("kicked- 
out" stars, jZolotov et al.||20T0l |Purcell et al.|[20T0l |Sheffield| 
|et al.|2012| ). Recent observations of halo stars (both local and 
distant field samples) provide evidence for the existence of 
at least two overlapping stel lar halos with distinct metallicity 
and rotation signatures (e. g. Carollo et al. 200 7) |2010| [Beers 
|eTaLl|20T2l |Hattori et aLl|2013[ |Kafle et al.||2012| ), however 
concerns have been raised about the d istances measurements 
used in some of these analyses (Schoenri ch et al.||2010| ) and 
their dependence on the assumed rota tion velocity of the local 
standard of rest (Deason et al. 2011). We do not discuss the 
possible formation modes of halo stars in this paper, but refer 
readers to other recent works listed above. 

The outline of this paper is as follows: in Section 2 we de- 
scribe the simulations used in this study and the methods used 
to analyze the orbits of stars and dark matter p articles drawn 
from the simulations (Laskar frequency analysis [La skar 1993; 
Vall uri & Merritt|[T998l |Valluri et al.||20T0| |2U12» . In Sec- 



tion |3.1| we measure the shape of the dark matter (and stellar) 
halo as a function of radius. We also analyze the types of halo 
orbits in both halos and characterize thei r sha pes and eccen- 
tricities as a function of radius. In Section |X2| we compare the 
phase- space distributions of halo stars and dar k-ma tter parti- 
cles using frequency maps. Finally in Section [33] we assess 
how baryonic condensation alters the orbits of the tidal debris 
that constitutes the stellar halo by examining the metallicity - 
formation time relation for halo stars and their orbital proper- 
ties. In Section|4]we summarize our results, and discuss their 
implications. 

2. SIMULATIONS AND NUMERICAL METHODS 

2.1. Simulations 

The MUGS sample of galaxies were produced using an N- 
body+SPH simulation of the formation of galaxies in a fully 
cosmological context. Dark matter halos of mass ~ 10 12 M 
at z = were identifie d in a dark-matter-on ly (TV-body) ACDM 
(WMAP3 parameters Sperg eTet al.|2007| ) simulation within a 
box of comoving length 50 /h Mpc, without regard to halo 
spin or mass accretion histor y. These halos were then re- 
simulated using GASOLINE (Wadsley et~aLl|2004| ) at high 
resolution using the volume renormalization technique to al- 
low for high resolution in the region of interest while simul- 
taneously including the full tidal field. The simulations in- 
clude metal cooling, heating from UV background radiation, 
star formation, stellar energetic and metal feedback, and metal 
diffusion. Star formation and stellar feedb ack are included via 
recipe s based on the "blastwave model" (Stin son et al.| 2006, 
2010). In the high resolution region of these simulations, a 



dark matter particle has a mass of 1.1 x 10 6 M , a gas par- 
ticle has an initial mass of 2.2 x 10 5 M , and star particles 
form with masses of < 5.5 x 10 4 M . All particles have a 
gravitational softening length of e = 312.5 pc. Full details of 
the MUGS simulations and properties of the sample galaxies 
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can be found in Stinso n et al. p010). The disk galaxy ana- 
lyzed in this paper (g 15784) is the largest of the disk galaxies 
from the MUGS sample. 

The galaxy g 15784 contains a bulge which is fitted by 
a de Vaucouleurs R 1 / 4 profile with an effective radius of 
r e = 0.67 kpc, and a disk with an exponential scale length 
of 3.38 kpc. The maximum halo circular velocity is 
360 km s -1 and it s virial radius (r vir ) is 29 0.8 kpc. We use the 
definition of r vir ( [Bryan & Norman|r9 98 ) which corresponds 
approximately to the radius within which the average density 
of the halo is ~ 100 times the critical density]^] This galaxy 
experienced its last major merger at z = 2.5 about 10 Gyr prior 
to the time at which we examine the simulation. 

Cosmological simulations predict that the dark matter halos 
of all galaxies contain hundreds to thousands of dark subha- 
los ( [Moore etal.|1998| ), although only about 10% of the total 
mass of the halo is in such subhalos. Subhalos in g!5784 
were identified by the AMIGA ha lo finder AHip1(Gill et al. 
[20041 |Knollmann & Knebe|[2009] ). To determine the frozen 
potential in which to integrate orbits of stars and dark mat- 
ter particles in g 157 84, we included mass within 4 virial radii 
to avoid an unphysical discontinuous density change beyond 
the virial radius. This region included several small satellite 
galaxies and numerous dark subhalos and a small mass con- 
tribution at the edge of the simulation volume from a neigh- 
boring 10 12 M galaxy ( [Nickerson etal.|201l|). Both the cen- 
tral stellar spheroid and luminous subhalos are unrealistically 
massive due to overcooling. Both dark and luminous satellites 
can perturb orbits of halo particles. Subhalos found by AHF 
were removed; however, a few subhalos of mass ~ 10 7 M 
remained in the inner regions of the galaxy, where automated 
halo finders have the m ost difficulty distinguishing substruc- 
ture ( [Knebe et al.|2011| ), making it difficult to fully assess the 
effects of the subhalos on the orbits of stars and dark mat- 
ter particles. We tested the characteristics of orbits integrated 
both with and without the gravitational effect of subhalos in- 
cluded in the frozen potential and found a higher fraction of 
chaotic orbits when subhalos were included. We also carried 
out co ntrolled simulations w here subhal os we re allowed to 
evolve ( [Debattista et al.|2008| — hereafter D08). We find both 
frozen and evolving subhalos increase the fraction of chaotic 
orbits. However the population that is most strongly affected 
is the resonant orbit family. These orbit which are "thin" (i.e. 
they occupy only a small volume of physical space), were 
shown to be important in halos where disks grow adiabati- 
cally ( |V12| ). When subhalos are "frozen" in place, a resonant 
orbit may repeatedly encounter a subhalo that lies close to that 
volume causing it to become chaotic. In this paper we only 
consider orbits integrated in the frozen potential from which 
subhalos more massive than 10 7 M were removed. This does 
not reduce the potential for scattering from smaller subhalos 
or irregularities such as tidal streams however, and we will see 
that this scattering probably contributes to the large fraction of 
chaotic orbits. 

We compare the orbits of particles from the cosmological 
simulation to orbits drawn from two controlled simulations 
of dis k galaxies gr own adiabatically inside dark matter halos 
(from |Vl^l & |vT2l ). 

In the controlled simulations the disks form quiescently and 
hence alter the halos adiabatically. They also do not have sub- 

7 Note that in VI 0, VI 2, r2oo was used. 

8 http://popia.ft.uam.es/AMIGA/ 



halos and other irregularities arising from hierarchical struc- 
ture formation. The first controlled simulation SA1 (see D08, 
|V10| & [VT2| for details) is a triaxial dark matter halo in which 
a disk of particles representing a baryonic mass fraction of 
2.5% was grown (starting from nearly zero initial mass) adia- 
batically and linearly on a timescale t g = 5 Gyr with the disk's 
symmetry axis parallel to the short- axis of the halo. The 
triaxial TV-body halo used in simulation SA1 was produced 
by multi-stage merger of spherical NFW halos. In the sec- 
ond controlled simulation, SBgs, a baryonic stellar disk forms 
self-consistently from 2.8 x 10 11 M of initially hot gas (con- 
stituting 10% of t he tot al mass) in a prolate NFW halo ( |V12| 
Debattista et al.| (|2013|)). In SBgs the hot halo gas initially 



has the same spatial distribution as the dark matter and slowly 
cools to form a disk with symmetry axis along the short axis 
of the halo. The halo and gas particles are given an initial 
specific angular momentum j, determined by overall cosmo- 
logical spin parameter A = (j /G)(|£|/M 3 ) 1/2 = 0.03 9, which 
is motivated by cosmological Af-body experiments (Bullock 
et al.||2001|). The simu lation closely follows that described 
in |Roskar et al.| ([2008 ) and is evolved with the par allel N- 
body+SPH code GASOLINE ( [Wadsley et aTl[2004| ) for 10 
Gyr. Table [T] summarizes some of the key properties of the 
different simulations used in this paper. 



2.2. Halo Orbit Sample 

In each simulation we selected 1 - 2 x 10 4 halo particles, 
randomly distributed within some volume. The "global sam- 
ple" refers to particles selected at random within a spherical 
radius r g centered on the model's galactic center (potential 
minimum). 

In the MUGS galaxy g 157 84 we identify halo stars by 
excluding stellar orbits associated with the massive central 
spheroid (orbits with apocenter radius r apo < 5 kpc) and stel- 
lar disk (orbits whose maximum distance from the disk plane 
Zmax < 3 kpc ). The global sample comprised ^ 10 4 halo stars 
with r g < 50 kpc. Since the controlled simulations do not form 
stellar halos, we follow the orbits of dark matter particles se- 
lected with r g < 200 kpc in these two simulations. 

In all three potentials orbits were integrated forward from 
their z = initial conditions for a duration of ~ 50 Gyr in a 
frozen potential resulting from the full mass distribution of 
the simulation including dark matter and baryons using an 
integration scheme based on the PKDGRAV tree. The as- 
sumption made here (and by most studies that involve the 
kinematics of halo stars in Milky Way) is that the potential 
of the galaxy at z = is nearly in steady state. While this as- 
sumption is not strictly true, there is evidence to support the 
view that halos of mass similar to the Milky Way, evolve less 
at the present epoch than they have in the past. |Busha et aL 
( |2007| ) show that halos with M 20 o < 10 13 M in ACDM sim- 
ulations evolved for 64 h~ l Gyr after the Big Bang (see their 
fig 8) grow rapidly till t/h = 10 (corresponding to 13.8 Gyr 
after the Big Bang) after which there is little change in their 
virial mass. This does not mean that secular evolution in the 
disk and halo stops at z = 0, merely that rate of mass accretion 
slows down. We emphasize that the long orbital integration 
times only serve to improve the accuracy with which orbits in 
the potential at z = can be characterized. We do not imply 
that the potential will have been static for the duration of inte- 
gration nor imply that the long integration time is physically 
meaningful. 
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Table 1 

Details of simulations 



Run 

Name 


r v ir 
[kpc] 


M vi F\ 

[10 12 m g 


] [10 11 M ] 


[pc] 


A0 

[kpc] 


Run Description 


Reference^] 


gl5784 

SAl 

SBgs 


290.8 
491.9 
317.9 


1.4 
6.9 
1.9 


2.1 

1.75 

0.90 


312.5 

100 

100 


3.38 

3.0 

1.9 


MUGS disk galaxy 

Triaxial halo+short-axis stellar disk 

Prolate halo+hot gas — > disk w/ SF 


[Stinson etal.|J2010) 

tto;[vT2] 

Debattista et al. (in prep) 



a M V i r : mass within the virial radius (r w y). 

mass in baryons. 
c Q)m: gravitational softening length 
d h: exponential radial scale length of disk. 
e Reference where the simulation is described. 
f Mass of baryons in stars only. 

2.3. Laskar Frequency Analysis 

In this section we provide a brief description of the Laskar' s 
frequency analysis method. A more detailed description of the 
application to orbits in Af-body simulations is given in |V10| 
and lVHl 

Orbits in galaxies ar e approximately quasi-periodic ( [Bin 



ney & Tremaine 2008; hereafter (BT08)), hence their space 



and velocity coordinates can be represented by a time series 

of the form: x(t) = X^=7 Ake luJk \ and similarly for y(t),v x (t) 
etc. This means that a Fourier transform of such a time se- 
ries will yield the spectr um of frequencies o^, and ampli tudes 
A k , that define the orbit ( jBinney & Spergelj 1982] [1984). In a 
three dimension potential, only three frequencies in the spec- 
trum are linearly independent. The remaining Uk are given by 
uj k = l k Q l -fmA + wA- For a regular orbit the frequencies 
Sli,i= 1...3 are constant "fundamental frequencies" that are 
related to 3 angle variables 6((t) - Vtjt and 3 integrals of mo - 
tion (more precisely actions J t ) (Binney & Tremaine 2008 ). 
Since the full phase- space trajectory in the action-angle co- 
ordinates (//,#;) resembles a 3D torus, the mapping from 
(x(f), v(Q) <^> (J, 6) is referred to as torus construction. Laskar 
( 1990| [T993] ) developed a very accurate numerical technique 
"Numerical Analysis of Fundamental Frequencies" to recover 
frequencies by taking Fourier Transforms of complex time se- 
ries of the form x(t) + i \(t). We use the implementation of 
Valluri & Merritt (1998), which uses integer programming to 
recover orbital fundamental frequencies with extremely high 
accuracy in ~ 20 - 30 orbital periods (significantly faster than 
100s of orbital period s needed by other codes; e.g., Carpin- 



tero & Aguilar 1998 Orbits were integrated and orbital fre 
quencies computed in Carte sian c oordinates and in cylindrical 
coordinates, as described in |V12[ 

Orbital fundamental frequenc ies can be used to: (i) classify 
orbits into major orbit families (Carpintero & Agu ilar|1998| ); 
(ii) quantify the elong ation of an orbit relative to the shape 



of the potential dVTOj); (iii) identify chaotic orbits (Laskar 
1993|); (iv) ident ify important resonant orbit families ( Robutel 
& Laskar 2001 ) which have trapped a large number of reso 



nant orbits yielding insights into the importance of secular 
or adiabatic processes; and (v) determine the self-consistent 
equilibrium potential from which a population of orbits was 
drawn ( |VT2] ). 

To determine the fraction of chaotic orbits in a potential, 
the orbital time series is divided into two equal segments and 
the orbital fundamental frequencies computed in each time 
segment. Since regular orbits have fixed frequencies that 
do not change in time, the change in the frequency mea- 
sured in the two time segments can be used to measure the 
drift in frequency space. |V10| showed that even for orbits 



in A/-body potentials (which are inherently noisy), it is pos- 
sible to distinguish between Af-body jitter and true chaos by 
a quantitative measurement of frequency drift by defining a 
diffusion rate parameter log(A/) which measures the loga- 
rithm of the change in the frequency of the leading term in 
the orbit's frequency spectrum, measured in two consecutive 
time segments. |V10| showed, using orbits in A/-body simu- 
lations of spherical NFW halos with gravitational softening 
e = 100 pc, that orbits with log(A/) < -1.2 were regular. In 
gl5784 the fraction of orbits with log(A/) > -1.2 (chaotic by 
the criterion established in the spherical NFW halo) is over 
50%. However visually many of the orbits have almost reg- 
ular appearances, g 157 84 has a larger gravitational softening 
(e = 312.5 pc) which appears to result in a slightly smoother 
potential and hence orbits with log(A/) < -0.5 appear more 
"regular" (i.e. they fill a well defined volume in configuration 
space) than they do in the controlled simulations. We there- 
fore use this more relaxed criterion for identifying regular or- 
bits in the cosmological simulation, but it must be kept in mind 
that in N-body simulations all orbits are inherently "noisy" 
resulting in a continuous distribution of diffusion parameters 
log(A/) and any cut-off between "regular" and "chaotic" 
is arbitrary. Furthermore we have found that decreasing the 
cutoff to log(A/) < -1.2 does not significantly affect the frac- 
tions of short- and long-axis tubes, but decreases the fraction 
of box orbits while increasing the fraction of chaotic orbits. 
We therefore caution readers that in this cosmological sim- 
ulation the boundary between chaotic and regular box orbits 
is rather fuzzy and although we distinguish between them for 
consistency with our previous work, they are better thought of 
as a single family of "centrophilic" orbits with no net angular 
momentum. 

3. RESULTS 

In this section we examine various properties of the stellar 
and dark matter halos of the simulated disk galaxy gl5784. In 
particular we e xamine the shape of the dark matter and stellar 
halos (Section [3.1.1 ) the orbital properties of the stars and 
dark matter particles (Section |3.1.2| ), and the phase spac e dis - 
tributions of halo stars and dark matter particles (Section |X2| ). 
It is important to keep in mind that although we expect that 
many of these results are likely to be typical of dark matter 
and stellar halos from cosmological hydrodynamical simula- 
tions, we are in fact, analyzing only one halo in a cosmologi- 
cal context and so the trends and conclusion found should be 
taken with caution since they might depend on the accretion 
history of the particular halo and on the sub-grid physics. 

3.1. Halo shape and orbital properties 
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Figure 1. Axis ratios b/a (top) and c/a (middle), triaxiality T (bottom) ver- 
sus radius from the center of galaxy g 15784 for: dark matter particles (solid 
lines) and halo stars (dot-dashed lines). The dark matter halo is oblate in the 
inner ~ 15 kpc transitioning to triaxial (0.33 < T < 0.66) at larger radii. The 
stellar halo is nearly oblate within within ~ 50 kpc (T ~ 0.2) and becomes 
slightly triaxial beyond this radius. 



3.1.1. Halo shape 

We computed the shapes of the spatial matter distributions 
of dark matter and star particles in g 157 84 within ellipsoidal 
shells centered on the minimum of the potential. The shape is 
determined at each radius r by finding the ellipsoidal shell 
of width ~ 5 kpc whose principal axes a, b, and c self- 
consistently (1) have a geometric mean radius y/abc = r, and 
(2) are the eigenvectors of the second moment tensor of the 
mass distribution within the ellipsoidal shell (J. Bailin et al. 
2013, in preparation). In practice, this is done by starting 
with a spherical shell, calculating the second moment tensor, 
deforming the shape of the shell to match the eigenvectors of 
the tensor, and iterating until the solution has converged; thi s 
is very similar to the method advocated by Zemp et al. ( 201 1\ . 

To measure the shape of the stellar halo (at z = 0) we select 
halo stars by excluding all stars within 5 kpc of the center of 
the galaxy (defined as the "bulge"), as well as all stars which 
lie within 3 kpc of the disk plane (defined as a thin+thick disk) 
for particles within 20 kpc (the radius within which most disk 
particles reside). This cut to identify "halo stars" is not perfect 
since it incorrectly excludes the small fraction of halo stars 
that happen to lie close to the disk plane or within 5 kpc of 
the galactic center at z = 0. This cut may potentially bias c/a 



values to be somewhat higher than they should be within the 
20 kpc inner region of the halo, but is unlikely to significantly 
alter b/a. We use a more sophisticated method (relying on 
kinematic properties) to exclude disk stars when we examine 
orbital properties of stars. 

Following standard nomenclature a,b,c are defined as the 
semi-axis lengths of long (x), intermediate (y) and short (z) 
axes respectively. The short axis is perpendicular to the plane 
of the disk in the inner region (the disk warps beyond about 
30 kpc). Figure [T] shows axis ratios as a function of radius b/a 
(top panel) and c/a (middle panel), and the triaxiality param- 
eter T = (l-b 2 /a 2 )/(l-c 2 /a 2 ) (bottom panel). The triaxiality 
parameter T is often used to characterize the shape of ellip- 
soidal figures which are termed "oblate" if T < 1/3, "triaxial" 
if 1/3 < T < 2/3 and "prolate" if T > 2/3. Thin dashed hor- 
izontal lines marking these shape regimes are marked in the 
bottom panel of Fig[l] 

For the dark matter halo, b/a and c/a vary only slightly 
with radius (1 < b/a < 0.9, c/a ~ 0.8). The bottom plot 
shows that the halo is oblate within ~ 15 kpc, becoming in- 
creasingly triaxial to prolate-triaxial beyond this radius. The 
shape of the stellar halo is highly oblate (T ~ 0.1) within 
40 kpc, becoming mildly triaxial at larger radii. The stellar 
halo is also much flatter (c/a is smaller) than the dark matter 
halo at all radii. Note that at radius ~ 17 kpc, the sharp dip in 
b/a and the corresponding spike in T are the result of a sub- 
halo of mas s < 10 7 M Q that w as not identified by the Amiga 
halo-finder. |Zemp et aL] ( |2011| ) show that when subhalos are 
not subtracted (or improperly subtracted) they produce spuri- 
ous spikes in the shape parameters. The overall changes in the 
shape parameters (e.g. T) of the dark matter halo and stellar 
halo as a function of radius are consi stent with the results of 
previous studies (e.g. Zemp et aL| 2012 ). 



3.1.2. Orbital type distributions 

To assess how orbital structure depends on shape, we ex- 
amined orbits of ~ 10 4 star and dark matter particles in the 
global sample of the galaxy g 15784, and for comparison, 10 4 
dark matter orbits from each of the two controlled simula- 
tions SA1 and SBgs. Orbits were classified into four types 
(box orbits, short-axis tubes [SAT], long axis tubes [LAT], 
and chaotic orbits). Their distributions as a function of orbital 
pericenter radius (r per i) and orbital eccentricity were examined 
in order to understand how the shapes of the two components 
(dark matter and halo stars) are related to the types of orbits 
that self-consistently reproduce their shapes. 

Figure [2] shows orbit populations as a function of r per i in 
the two controlled simulations SA1 and SBgs (top panels), 
and in g 157 84 (bottom panels). We choose to plot the orbit 
fractions as a function of r per i rather than some other measure 
of the "radius" of an orbit (such as the time-averaged radius) 
because we are specifically interested in understanding if the 
SATs contribute to the oblateness in the inner halo. Since r per i 
defines the inner boundary of a tube orbit, it defines the region 
interior to which the orbit contributes no mass. Hence this 
quantity is best suited to addressing this question. In theory, 
box orbits and chaotic orbits can travel arbitrarily close to the 
center of the potential and so formally they have r per i = but, 
in practice r per i > during the 50 Gyr (and at least 20 orbital 
period) integration time. Each panel shows kernel density dis- 
tributions of each of the four orbit families as indicated by the 
line legends. The area under each curve is proportional to the 
total number of orbits of that family. It is important to note 
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Figure 2. Kernel density distributions of orbit types as a function of orbital pericenter distance r per i for 10 4 halo particles in 3 simulations. In each panel the 
vertical axis is proportional to the number of orbits with a particular value of r per i, i.e. the sum of integrals over all curves in a panel is equal to the total number of 
orbits. Top left: orbits of dark matter particles from adiabatic A/-body simulation SA1 (triaxial dark matter halo with stellar disk). Top right: orbits of dark matter 
particles from adiabatic simulation SBgs (prolate dark matter halo with stellar disk formed from hot gas). Bottom left: dark matter particles in cosmological 
simulation g 15784. Bottom right: halos stars in g 157 84. Line styles indicate the 4 major orbit families. Only orbits with periods shorter than 2.5 Gyrs were 
selected. In the top two plots the spatial gravitational softening length is 0.103 kpc, while in the bottom two panels it is 0.3125 kpc. The orbital distributions 
below the softening radius in each model are unreliable. 



that in the top two panels the gravitational softening length of 
dark matter particles is ~ 0. 1 kpc, while in the bottom panels 
the gravitational softening length is ~ 0.3 kpc; distributions 
at pericenter radii smaller than the softening lengths are not 
reliable. 

The progenitor halo of SA1 was triaxial with 86% of parti- 
cles on box orbits and 11% on LATs dV10| ). Figure [2] shows 
that this population was transformed (due to the growth a stel- 
lar disk) to a halo that is still dominated by boxes (49%) but 
the fraction of SATs (32%), LATs (12%), and chaotic orbits 
(7%) has increased significantly. In contrast the progenitor 
halo o f SBgs was prolate with 78% LATs and 15% box orbits 
( |V10| ). Following the condensation of hot gas and the growth 
of a stellar disk, the dark matter halo orbit population evolves 
to an overall potential in which SATs dominate (51%) with 
LATs (35%) also important. In the cosmological simulation 
g 15784, halo star particles (Fig. [2] bottom right) are primar- 
ily on SATs (47%) and LATs (29%) and dark matter particles 
also have similar orbit fractions (40% are on SATs and 36% 
are on LATs). 

Table [2] shows the orbit fractions in the inner and outer re- 
gions of all three halos following the growth of disks. We 
computed the relative fractions of different types of orbits 
with r per i less than or greater than 20 kpc. Since the inner 
~ 20 kpc region is where the disk dominates, this is the re- 
gion where the shape of the halo is close to oblate axisym- 
metric in all three potentials. In none of the halos do SATs 
constitute more than 50% of the orbit fraction in the oblate in- 
ner halo (r per i < 20 kpc). In all four models the "triaxial orbit 
families" (box orbits and chaotic orbits and long-axis tubes), 



together constitute 50% or more of the population. In SBgs, 
g 15784, the outer part of the halo (r per i > 20 kpc) are mildly 
triaxial. Here we find that LATs are the dominant popula- 
tion. In SA1 - which is quite strongly triaxial beyond 20 kpc 
- SATs constitute 65% of the population with LATs constitut- 
ing 33%. Thus in no case do we find evidence that an oblate 
density distribution implies a distribution function dominated 
by axisymmetric SATs, and conversely we find that the pres- 
ence of a large fraction of short-axis tubes does not imply 
axisymmetry. 

Although the distributions of orbit types as a function of 
r per i differ in detail for dark matter particles and halo stars in 
g 15784, their overall distributions are significantly more sim- 
ilar to each other than they are, for instance, to the orbit distri- 
bution of dark matter particles in model SA1 (Fig.[2]top left). 
Furthermore g 15784 shows a striking resemblance to the or- 
bit populations in SBgs. The similarity is probably coinciden- 
tal in view of the fact that g 15784 formed in a cosmological 
context and has experienced significant perturbations due to 
mergers and substructure. Minor mergers are entirely absent 
from the history of SBgs, which formed from a single major 
merger between two spherical halos. We expect the halos of 
galaxies of the mass of g 157 84 to have experienced multiple 
major mergers and to therefore be more like SA1, rather than 
SBgs which has only had one major merger. It appears that 
since gl5784 had its last major merger nearly 10 Gyr ago, this 
galaxy's stellar halo was primarily formed via minor majors 
(rather than major mergers). Since the orbits at large radii are 
primarily tubes, we hypothesize that the satellites that formed 
the stellar halo were mainly accreted on tube orbits. We will 
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Table 2 

Orbit fractions in inner and outer halo regions 



Run 




r p eri < 20kpC 






r p eri > 20kpC 






Box 


Chaotic 


LAT 


SAT 


Box 


Chaotic 


LAT 


SAT 


SA1 


0.51 


0.09 


0.11 


0.29 


0.01 


0.05 


0.33 


0.65 


SBgs 


0.10 


0.10 


0.30 


0.50 


0.03 


0.02 


0.52 


0.43 


Stars 


0.12 


0.13 


0.32 


0.43 


0.03 


0.02 


0.55 


0.41 


DM 


0.12 


0.12 


0.34 


0.42 


0.09 


0.00 


0.59 


0.32 



see in Section|3.3| that there is su pport for this hypothesis. 

Recently Brya n et aL] ( |2012| analyzed orbits in halos of 
mean mass ~6x 10 1J M© (at z = 0) and 7 x 10 11 M (at z = 2) 
from the Overwhelmingly Large Simulations (OWLS) to in- 
vestigate the effects of various feedback prescriptions on the 
orbital properties of stars and dark m atter particles in a cos- 
mological context. Bryan et al. (2012) selected 50 halos from 
5 different simulations (with different feedback prescriptions) 
and selected 500 particles per halo. They integrated orbits 
in gravitational potentials computed with the Self-Consistent 
Field (SCF) method and cla ssified them using the method of 
Carpi ntero & Aguilar|fl998| ). 

|Bryan et al.| ( |2012f find that as the central baryonic frac- 
tion increases, the fraction of box orbits in the inner part of 
the halo decreases. Since we only examine one simulation 
we are unable to verify their result that the fraction of box 
orbits in the inner region decreases with increasing central 
baryonic mass fraction; however, the baryonic mass fraction 
in the inner 20 kpc region of g 15784 (where the disk domi- 
nates) is obviously much higher than in the outer region, and 
we find no evidence that box orbits have been depleted in 
the inner regions as a consequence of baryonic c ondensation 
The d ifference between our result and those of |Bryan et al. 
( |2012| ) is most probably a consequence of differences in the 
way in which the total galactic potential (in whi ch orbits are 
evolved) is computed. The SCF method (used by Bryan et al. 
2012| ) uses a low-order multipole expansion code to compute 
the potential, which consequently is much smoother and has 
more symmetry than our potential, which is computed with 
the PKDGRAV tree-code and therefore includes all the irreg- 
ularities of the full cosmological simulation (except subhalos 
more massive than 10 7 M which were removed). 

3.1.3. Orbital shapes: elongation versus eccentricity 

In any self-consistent potential, the shapes of the major- 
ity of the orbits must match the overall shape of the 3- 
dimensional matter density distribution. In a triaxial poten- 
tial, the elongation along the major axis is provided either by 
box orbits or inner LATs. The ratios of the fundamental fre- 
quencies of orbits can be used to characterize their overall 
shape. For a triaxial potential, the fact that semi-axes lengths 
a > b > c implies that the oscillation frequencies obey the 
condition | O x | < | Q y | < \Q Z \ for any (non resonant) orbit with 
the same over-all shape as the density distribution (we con- 
sider only the absolute values of the freq uenci es since their 
signs only signify the sense of oscillation). |V10| use this prop- 
erty to define an average "elongation" parameter Xs for any 
orbit. When an orbit is elongated in the same way as the po- 
tential, 

\n z \ > \n y \ > \n x \ 

|«z| l«zl" 



Therefore we can define an elongation parameter Xs by 



such that Xs > for orbits elongated along the major axis with 
larger values of Xs implying a greater the degree of elonga- 
tion (V10). |V12| showed that some systems can be elongated 
along the y-axis at small radii but elongated along the x-axis at 
larger radii; also some outer LAT orbits can have greater ex- 
tent along the y-axis than along the x-axis. In both these cases 
Xs is slightly negative. Orbits for which all frequencies are 
almost equal enclose a volume that is almost spherical have 
Xs ~ because Q x ~Q y ~Q z , Also orbits that are close to ax- 
isymmetric about the short (z) axis (i.e. SATs) have Q x ~ Sl y 
and hence Xs ~ 0, regardless of the value of Q z . 

Figure [3] shows kernel density distribution of orbits of the 
four different types as a function of orbital elongation Xs- 
SA1 shows significantly larger fraction of elongated orbits 
(Xs > 0.1) than other models. The SATs (dashed lines) have 
Xs = because they are close to axisymmetric. Model SBgs 
(which is prolate at large radii) has a significant fraction of 
elongated LAT orbits. In contrast, for g 15784, all four fam- 
ilies of orbits for both dark matter (bottom left) and stars 
(bottom right) have very small elongation (i.e. they are quite 
"round" or axisymmetric). 

|V10| showed for several controlled simulations (including 
SA1) that the orbital elongation parameter Xs in the inner re- 
gions decreased from Xs ~ 0.4-0.5 in the triaxial or prolate 
models to Xs ~ 0-2 after baryonic components were grown. 
They showed that although some orbits do change from boxes 
or LATs in the prolate systems to SATs, all orbits adapt to 
the more axisymmetric baryonic component by becoming less 
elongated, especially at small r per i. However, only a small 
fraction of all orbits (between 4% and 25%) actually trans- 
form themselves to become SATs. Figure [4] shows the median 
value of Xs in 15 bins in r per i for orbits m the three simula- 
tions. The top two panels show that when simulations are 
controlled (adiabatic) the baryons make box and chaotic or- 
bits with smaller pericenter radii less elongated but at larger 
pericenter radii these orbits remain elongated (LATs always 
tend to be elongated but as shown by |V10 , these too are more 
elongated when baryons are absent). The lower two panels 
show the elongations of orbits in g 157 84. In this simulation 
all types of orbits are essentially "round" over the en tire range 
in r per i. Although this confirms the findings of |V10| that orbits 
of all three families can become "round" to support the shape 
of the potential, it is also clear that overall the shapes of all 
the orbits are similar and close to axisymmetric at all radii. 

Another commonly used metric for characterizing individ- 
ual orbits is the eccentricity parameter 

^ _ ^apo — ^peri 
F apo ~l~ 7" peri 

Strictly speaking, orbital eccentricity e is defined for a two 
dimensional orbit in terms of its apocenter radius r apo and 
pericenter radius r per i. It is important to emphasize that ec- 
centricity does not measure the degree of axisymmetry; e.g. 
a perfectly planar orbit in a spherical potential is always an 
axisymmetric rosette, but can have an arbitrary eccentricity 
that depends on its angular momentum. At a given energy, 
orbits with very high angular momentum have eccentricity 
approaches zero while orbits with very low angular momen- 
tum have eccentricity approaches unity. Figure [5] shows the 
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kernel density distributions of orbital eccentricity for parti- 
cles in the three simulations. We see that all four orbit fami- 
lies have larger fractions of orbits with high eccentricity than 
with low eccentricity. Even SATs (dashed curves) and LATs 
(dot-dashed curves) are very eccentric implying that they have 
fairly low angular momentum. The absence of low eccentric- 
ity box and chaotic orbits is not unexpected - these orbits are 
highly radial (formally having zero time averaged angular mo- 
mentum) and have r per i approaches zero. The average orbital 
eccentricities of halo stars in each orbit family are e = 0.88 
(boxes), e = 0.95 (chaotic), e = 0.76 (SATs) and e = 0.72 
(LATs). For dark matter particles the average orbital eccen- 
tricities of each orbit family are e = 0.82 (boxes), e = 0.95 
(chaotic), e = 0.71 (SATs) and e = 0.63 (LATs). These average 
orbital eccentricities for the various orbit families are some- 
what larger than the average orbital eccentricities (e = 0.6) of 
dark matt er subhalos and particles in dark-matter-only simu- 
lations ( van den B osch et al.|1999| ). It is particularly interest- 
ing to reflect on the fact that for g 15784 halo stars, even those 
on SATs, are more likely to have high eccentricity than low 
eccentricity. If a significant fraction of the stellar halo was 
formed in situ in an early disk that was then heated to form 
the halo, one would expect this population to have a higher 
angular momentum and lower eccentricity than an accreted 
population. The fact that the majority of the stellar halo or- 
bits in this global sample have a relatively high eccentricity 
reflects the fact that they were primarily accreted. 
The main conclusions from Figures [2] -[5] are: 

• The types of orbits that characterize a halo (either in a 
controlled simulation or in a cosmological simulation) 
reflect both the merger/formation history and the effect 
of baryons. 

• All four types of orbits (boxes, SATs, LATs and 
chaotic) can be present in significant proportions even 
when the shape of the potential is nearly oblate. 

• All four types of orbits are able to become less elon- 
gated (i.e. more axisymmetric) to adapt to the nearly 
oblate shape 

• All four types of orbits have a high fraction of eccentric 
orbits reflecting the fact that the halos are products of 
mergers or radial infall. 

3.2. Phase space distributions: Halo stars versus dark matter 

Since the mass of the stellar halo is only a tiny fraction of 
the mass of the dark matter halo, the dynamics of halo stars 
are determined by the shape of the dark matter halo, the orbital 
initial conditions at the time of accretion, and th e evolution of 
the ha lo potential since the epoch of accretion ( [Knebe et al.| 
|2005| ). All of these factors also influence the orbital types, or- 
bital shapes and eccentricity distributions of orbits. The sim- 
ilarities in the distributions of orbits of dark matter and halo 
stars with r per i, \s and eccentricity (Figs. [2] -[5} show quali- 
tative similarities that are quite significant, especially in com- 
parison with the orbit populations of the controlle d simulation 
SA1, which formed from multiple major mergers. Brya n et aL] 
(2012) also found in their analysis of orbits from the OWLS 
cosmological simulations that star particles and dark matter 
particles have similar orbital distributions as a function of ra- 
di us. 

|V12| showed that a useful way to map the entire distribu- 
tion function of a self-consistent system was via frequency 



maps. A frequency map is obtained by plotting the ratios 
of fundamental frequencies: in Cartesian coordinates Q x /Q z 
vs. Q, y /Q, z , or in cylindrical coordinates, vs. 0,^/0,^ 

They showed that frequency maps give a pictorial represen- 
tation of the different orbit families, their relative importance 
and how they are distributed with energy. In the frequency 
maps in Figure |6| each point represents a single orbit with 
color representing the orbital binding energies in 3 bins (blue: 
the l/3rd most tightly bound; red: l/3rd least tightly bound; 
green: intermediate energy). 

Figure [6] shows the frequency maps for orbits in galaxy 
g 15784. The top left panel shows a map (in Cartesian coor- 
dinates) for dark matter particles, while the bottom left shows 
the same for halo stars. Short-axis tubes lie along the di- 
agonal £t x /£t z ~ Qy/Q z - The long-axis tube family is clus- 
tered along the horizontal line Q, y /Q, z ~ 1, box and chaotic 
orbits are mostly the blue points scattered around the map. 
The dashed lines mark the major resonances (satisfying con- 
ditions like lQ x + m£l y + nQ z = 0), with numbers representing 
the integers l,m,n. In this Cartesian representation no orbits 
are seen to be associated with box-orbit resonances (diagonal 
lines on the left-side of the graph). The right hand panels of 
Figure [6] shows the frequency maps in cylindrical coordinates 
for the same sets of orbits. A prominent resonance is seen 
diiVL z /Q J R = Q J( f ) /VL R ^ Q.5 i n orbits of both the halo stars and 
dark matter particles. |V12| showed that numerous resonances 
characterize the phase space of a self-consistent distribution 
function evolved adiabatically in controlled simulations. This 
resonance is associated with a thin shell resonance of the SAT 
family and it is impressive that such resonant orbits persist de- 
spite the hierarchical accretion history of this system. The fre- 
quency map representations of the halo stars and dark matter 
particles are very similar to each other both in Cartesian maps 
(left) and in the cylindrical maps (right) providing additional 
evidence that their phase space distribution functions are sim- 
ilar. This implies that the phase space distribution function 
of this sample of stellar halo orbits probably originated in a 
manner similar to that of the dark matter orbits. 



3.3. 



Galactic archeology: age -metallicity- orbital 
correlations 



Since the prescient work of Eggen et al. 



fT962] ), it has been 

recognized that correlations between orbitalprbperties of halo 
stars and their metallcities and ages can be used to infer the 
formation history of the Galaxy. If the orbital integrals of 
motion of halo stars (e.g. total energy E, the total angular 
momentum L, and angular momentum about the z-axis L z ) 
change only slowly as the potential evolves and the debris 
disperses, correlations in integral-of-motion space can lead 
to the recovery of the progenitor satellite and also enable re- 
covery of the underlying potent ial from tidal debris even in 
time-dependent potentia ls ( Johns ton et al.|199 6; Helmi & de 
Zeeuw 2000; Bullock & Johnston 2005 ). Searching for cor- 
relations m angle-action or frequency space may improve the 
accuracy of such recovery ( McMillan & Binney 2008 ; Gomez 



|& He lmi 2010). Using semi-analytic prescriptions to con 
struct stellar por 



populations for accreted hal os from cosmolog 
ically motivated pure Af-body simulations |Font et aL] ([2006 ) 
showed that combining phase space information with chem- 
ical abundances and stellar ages can aid the recovery of in- 
dividual satellites, even when the stars in a satellite have a 
small spread in ages and chemical characteristics. To date 
there have been no studies which have tried to identify halo 
stars associated with individual satellite progenitors in a cos- 
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Figure 5. Similar to Fig [2] Kernel density distribution of orbit types as a function of their eccentricity. 



mological hydrodynamical simulation. 

Figure [7] shows orbital energy versus one component of 
angular momentum for halo stars from the globa l samp le in 
each of the four orbit families studied in Section 13.1.21 For 
box orbits, SATs and chaotic orbits the abscissa is the an- 
gular momentum about the z axis (L z ), while for LATs the 
abscissa is L x since x is the axis about which this tube fam- 
ily rotates. Points corresponding to orbits in each family are 
fairly smoothly distributed and almost no substructure typi- 
cal of tidal debris associated wi th distinct individual satellites 
(e.g. |Helmi & de Zeeuw||2000 ) is seen. The ~ 10 4 particles 
plotted in Figure [7| were selected at random within a spherical 
volume of r g = 50 kpc. However, even a local sample selected 
within a small volume around the location of a fiducial "sun" 
showed no substructure in £",L,L Z ,L X space. It is important to 
emphasize that we only plot E versus L for comparison with 
previous work but that in fact none of the components of angu- 
lar momentum are expected to be integrals of motion for any 
of the orbit families in this triaxial potential. Figure [7] clearly 
shows that with in this simulation, the hierarchical growth of 
the halo, accompanied by dissipative evolution, have erased 
all correlations between these pseudo integrals of motion. It 
should be noted that our choice of orbital initial conditions 
(randomly chosen within 50 kpc of the galaxy's center) is 
not ideal for searching for coherent tidal debris. Typically 
for such studies a local volume representing the solar neigh- 
borhood is examined. Our examination of such a volume (i.e. 
R s = 4 kpc yielded too small fraction of bona fide halo stars 
(because the disk in this simulation is quite thick) to perform 
this exercise. We also searched for substructure in plots of 
energy vs. metallicity, energy vs. stellar age and angular mo- 
mentum vs. metallicity and stellar age. None of these plots 
revealed any detectable substructure that would enable us to 
identify satellite progenitors. 

Bailin et al. (2013, in preparation) analyze the age- 



metallicity relation for all stars (halo and disk) in the MUGS 
cosmological disk galaxy simulation g 157 84. They find that 
while the most metal-poor stars in the galaxy were formed at 
the earliest times, the ISM is rapidly enriched and hence the 
formation of highly metal enriched stars follows very quickly 
thereafter. Bailin et al. show that when [Fe/H] is plotted 
versus formation time, numerous tight "stream-like" features 
which correspond to continuous star formation within individ- 
ual galactic subunits, are seen. A study of the relative spatial, 
metallicity and age distributions of stars that were accreted 
onto g 157 84 and stars that formed in the disk and were sub- 
sequently kicked into the halo will be presented in Bailin et 
al. Here we focus only on 10 4 halo stars whose orbits were 
analyzed in previous sections, and examine how their orbital 
characteristics depend on their metallicities and ages. 

In Figure [8] we plot total metallicity (in solar units) for stars 
in the halo orbit sample studied in previous sections, versus 
their formation time in the simulation (tf OTm ). The formation 
time is scaled such that stars that form at the present time 
(z = 0) have tf orm = 13.7 Gyr. The panels shows stars on differ- 
ent types of orbits: boxes, SATs, LATs and chaotic orbits as 
indicated by the labels. Each point represents a star particle 
from the simulation. 

Most striking are the thin stream-like features (identified 
by Bailin et al.) which arise from progressive enrichment 
and star-formation inside individual sub-galactic units. Note 
that even though we have randomly selected ~ 10 4 halo stars 
within 50 kpc, with no other constraints other that they have 
r apo > 5 kpc and Z max > 3 kpc the thin streams associated with 
individual progenitors are clearly visible. Star formation ap- 
pears to stop at some point in all the streams. While more 
detailed analysis with all the MUGS snapshots are needed to 
determine the precise reason for this, a plausible interpretation 
is that star formation in a subhalo can continue for some time 
after it is accreted by the galaxy but stops once it is tidally 
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Figure 6. Frequency maps for dark matter particles (top row), halo stars (bottom row) selected within 50 kpc of the galactic center. Left: Frequencies computed 
in Cartesian coordinates, Right: frequencies computed in cylindrical coordinates. Colors represent binding energies of particles and dashed lines mark major 
resonances (see text). 



disrupted. Consequently the end of the "stream" serves as 
a proxy for when it is significantly disrupted by tidal forces 
from the main galaxy. 

The thin stream-like features which were disrupted early 
(/form < 8 - 10 Gyr) appear in all four panels indicating that 
at z = 0, stars associated with a single progenitor satellite tra- 
verse the halo on all four types of orbits. For example the 
features labeled "A" consists of stars with a wide range of 
metallicity and tf 0Ym and are clearly seen in four panels. This 
implies that these two galactic sub-units were completely dis- 
rupted by tidal forces and then scattered onto different types 
of orbits by mixing processes. 

We believe this is probably "chaotic mixing" rather than 
"phase mixing" since the latter conserves orbital integrals of 
motion in a static potential and changes them adiabatically in 
a slowly varying one. Phase mixing is not capable of chang- 
ing th e orbital type, but chaotic mixing is (Merritt & Valluri 
|1996| ). Therefore it is unlikely that the phase space coordi- 
nates of stars (at z = 0) in the two units labeled "A" can be used 
to fully reconstruct the progenitor or its orbital properties at 
the time it was accreted. For dynamically older streams like 
"A" the integrals of motion, especially energy, are expected 
to change significantly due to the change in the halo potential 
due to mass accretion, both dark ( Knebe et al. 2005 ) and bary- 
onic (jVlO]). Additionally chaotic mixing in the triaxial halo 



causes tida l streams to disperse more rapidly ( Vogelsberger 
|et al.|2008D . 

In contrast the tf orm -metallicity feature labeled "B" is most 
clearly seen in the top right panel (SATs). Its stars also span 
a range of metallicities (—1.5 < [Z/Z ] < -0.5) and forma- 
tion epochs (6-13.5 Gyr) however most of the stars of this 
satellite are on SATs (although a faint hint of it is seen in the 
LATs panel) while there is no distinct counterpart in the box 
or chaotic orbit panel^] This indicates that the orbits of this 
satellite are not mixed, and hence the progenitor's phase space 
coordinates and the halo potential can probably be extracted 
by backward evolution of stellar orbits. The stream labeled 
"C" in the bottom left panel consists of stars on LATs and 
is also only seen in this panel. The near-exclusivity in orbit- 
type for streams "B" and "C" implies that the orbits of these 
stars have not experienced much mixing (probably because 
they were accreted/tidally disrupted more recently). Figure [8] 
shows that the thin steam-like features have stars with a wide 
range of metallicities and star formation occurs over an ex- 
tended period of time. We see that stars with a range of metal- 
licities and ages that were once associated with the same pro- 
genitor are easy to identify because they form thin features 
in metallicity-formation time plots despite sometimes having 
very different orbital characteristics. 

9 although admittedly this could be in part due to inadequate sampling. 
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Figure 7. Total energy versus one component of angular momentum for the 
global sample of halo stars in g 15784. Panels show stars on different types of 
orbits as indicated by the labels. For Box, SAT and chaotic orbits ordinate is 
L z , for LATs (bottom left) the ordinate is L x . 

Our preliminary analysis of the age-metallicity signatures 
of our selected halo stars in Figure [8] shows that satellites 
leave the stars in long and short-axis tube orbits as they ac- 
crete. Satellites that accreted >5 Gyr ago have chaotically 
mixed these stars into chaotic and box orbits, but more re- 
cently accreted satellites have not yet become well mixed. 
Such chaotic mixing of early accreted satel lites is consistent 
with previous studies of satellite accretion ([Bullock & John 
|ston|2005l|Font et al.|2006[|G6mez & Helmi|2010| ). " 

It is important to note that star formation appears to con- 
tinue in this halo population till quite late times — well be- 
yond what is expected based on observations of the stellar ha- 
los of the Milky Way and M3 1 . This is in part due an artifact 
of the detailed sub-grid physics in the simulations: the cool- 
ing of gas in to dark matter subhalos results in their becom- 
ing more massive than in the real Universe. This also allows 
subhalos in these simulations to continue to form stars after 
they are accreted by the main galaxy. This (presumed) arti- 
fact is not unique to our simulation s and has been previ ously 
observed ( |Tissera et al.|2012[|20T3] ). |Tissera et al.l ( [2013] ) refer 
to this component of the stellar halo as "endo-debris", and dis- 
tinguish it from "debris" stars which were formed in a satellite 
before it was accreted. There is tentative evidence that such a 
population has also b een discovered in the Milky Way halo by 



Sheff ield et al.| ( |2012] ) who have identified a population of halo 
stars with low [Fe/H] but high [a/Fe], which they argue could 
have been formed in si tu in the halo (i.e. inside satellites af- 
ter they were accreted). Tiss era et al. ( |2013| ) also suggest that 
on e population of Carb on Enhanced Metal Poor stars found 
by |Carollo et al.| ([2012) could be such "endo-debris" stars. If 
star formation in satellites terminates earlier than in these sim- 
ulations, narrow features would not extend to late times, but 
would still remain detectable at early times as in Figure [8] 



Together Figures [7] and [8] imply that the tidal debris in 
cosmological hydrodynamical simulations experiences signif- 
icantly more chaotic evolution than in collisionless simula- 
tions, making it much harder to identify individual progen- 
itors using phase space coordinates alone. But the identifi- 
cation of distinct progenitors is quite easy in age-metallicity 
plots, and such plots when combined with orbital information 
can be used to probe the accretion history of the halo. Cur- 
rent cosmological hydrodynamic simulations, including the 
one studied here, do not have adequate resolution to iden- 
tify the types of stellar streams observed in the Milky Way 
and M31. The study of the correlations between kinemat- 
ics, metallicities, and spatial distributions of halo stars and 
their formation history will become increasingly important in 
the future, as the resolution of simulations and the baryonic 
physics prescriptions improve. 

4. SUMMARY AND CONCLUSIONS 

We have analyzed the orbital properties of halo stars and 
dark matter particles from a disk galaxy forming in a cosmo- 
logical hydrodynamical simulation from the MUGS project. 
The simulated disk galaxy we have studied is a reasonably 
good Milky Way analog, although it has a central bulge and 
satellite galaxies that are too massive. In this galaxy the 
shape of the dark matter and stellar halos vary with radius due 
to the central condensation of baryons: oblate within about 
20 kpc, and mildly triaxial at i ntermediate and large radii, 
consistent with previous work ( Dubinski & Carlberg 1991; 
Kazantzid is et al.||2004[ JDebattista et al.||20081 |Tissera et al.| 
|2010[|Kazantzidis 'et al |2010||Zemp et al.|2012) . The stellar 
halo is slightly more oblate than the dark matter halo at all 
radii. We find that although the inner halo region - where the 
disk dominates - is oblate, it is dominated by box orbits and 
chaotic orbits and less than 30% of orbits in the inner oblate 
region are SATs. This is contrary to the general expectation 
that as an ellipsoidal figure becomes more oblate, i ts distri- 
bution function is increasingly dominated by SATs ([Ge rhard 
& Bi nney|1985T|Merritt & Valluri|1996HBinney & Trem aine 
2008). In the outer regions, where the potential is more elon- 
gated, SATs and LATs appear in roughly equal proportions 
and box orbits and chaotic orbits are insignificant. Thus in no 
case do we find evidence that oblate regions of the potentials 
are supported mainly by SATs as is generally assumed. Our 
findings confirm our previous work with controlled simula- 
tions (D08,|VT0l). 

The phase space distributions of halo stars and dark matter 
particles (as characterized by orbit populations as a function 
of radius and eccentricity) in the cosmological hydrodynami- 
cal simulations are very similar to each other. Although orbits 
of both types have a wide range of orbital eccentricities, the 
majority (even of those on SATs) have high eccentricity, in- 
dicative of their largely accreted origin; the distributions of 
orbital types with pericenter radius, orbital elongation param- 
eter and orbital eccentricity differ only slightly. Orbits of all 
four families have three dimensional shapes (as measured by 
the elongation parameter \s) that are close to axisymmetric, 
enabling them to self-consistently generate the nearly oblate 
shape of the inner halo. Where there are differences between 
the orbits of star and dark matter particles, they are clearly at- 
tributable to the differences in their spatial distributions. For 
instance, since the dark matter distribution is more extended 
and more triaxial than the stellar distribution, it is dominated 
by LATs in the outer regions. The less extended stellar halo 
is also more oblate and hence SATs dominate even beyond 
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20 kpc. Similarities in the orbital phase space properties (and 
distribution with energy) are also seen in the frequency maps. 

The similarities between the orbital structures of the stellar 
halo and dark matter halo suggest that both components share 
a common dynamical origin, probably accretion, although it 
is possible that some of the similarities arise because chaotic 
scattering of orbits of both types of particles, in the fully cos- 
mological simulation. The relatively large fraction of chaotic 
orbits may not be representative of real galaxies and is most 
likely a consequence of the presence of the massive central 
bulge and subhalos which produces some scattering (D08), 
and the "freezing" of irregularities like small subhalos and 
tidal features that accentuate the scattering of resonant orbits. 

Several previous studies have focused on uncovering tidal 
streams from the orbital phase space coordinates of halo stars 
in dissipationless Af-body simulations. In this fully cosmolog- 
ical hydrodynamical simulation we find that energy and an- 
gular momentum (or other quantities resembling integrals of 
motion such as orbital actions) — the variables traditionally 
used to identify tidal de bris in idealized collisionle ss simula- 
tions of tidal disruption ( He lmi & de Zeeuw|20 00) — evolve 
so significantly that they can no longer be used (exclusively) 
for identifying substructure in phase space. The importance 
of chaotic mixing in such simulations is evidently larger than 
in the controlled N-body simulations studied previously. 

We show that in plots of metallicity versus formation time 
(or stellar age), thin stream like features can be associated 
with individual galactic satellites in which star formation oc- 
curs continuously over a long period of time before it is tidally 
disrupted in the main galaxy's halo (Fig. [81. Consequently 
stars associated with a single disrupted satellite are neither of 
single metallicity nor do they have a narrow range of stellar 



ages (as also found by |Font et a l. 2006 ). For older disrup- 
tion events matched stream-like features are found on all four 
families of orbits implying that they have experienced chaotic 
mixing. However, younger features are typically associated 
with a single orbit family. 

The results of this paper demonstrate the value of orbital 
analysis for probing dynamical structure and formation his- 
tory. Similar studies of simulated triaxial galaxies and dark 
matter halos show that systems that form via dissipation- 
less " dry mergers" contain primarily bo x orbits and long-axis 
tubes ( [Barnes 1 1996H Valluri et al.|2010| ) while tube o rbits tend 
to be prod uced by dissipation in "wet mergers" (Hoffman 
|et al.|2010| ). Conversely if a system is dominated by box orbits 
or chaotic orbits one can infer that dissipationless formation 
such as a multi-stage merger or that scattering by a dense cen- 
ter has been important in the evolution of the system. 

It is important to point out that the MUGS simulations used 
in this paper have recently been superseded by several simu- 
lations in which new feedback prescriptions are able to pro- 



disk/total ratios ( 


Governato et al.||20101 Guedes et al.||2011l 


Zemp et al.|2012 


Stinson et al.|2012). Future generations of 



simulations will continue to improve the match between ob- 
servations and simulations. It is important to carry out similar 
analyses in these improved simulations to assess how specific 
implementations of baryonic physics alter the star formation 
rates and enrichment of subhalos that make up the stellar halo. 

In coming years phase space coordinates, elemental abun- 
dances and accurate isochronal ages for many halo stars will 
be obtained with Gaia and future wide-field, high-precision 
synoptic surveys, such as LSST. Surveys such as Gaia have 
been designed to yield data with observational errors of 10% 
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or less for stars within 10 kpc of the sun. Several ground based 
surveys will go even deeper. In a future paper we will con- 
sider the effects of selection bias and observational errors on 
the phase space coordinates, stellar ages, and elemental abun- 
dances on the recovery of the orbital structure of the stellar 
halo. 
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